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Abstract 

We present a new determination of the strong coupling constant from 
lattice QCD simulations. We use four different short-distance quantities 
to obtain the coupling, three different (infrared) meson splittings to tune 
the simulation parameters, and a wide range of lattice spacings, quark 
masses, and lattice volumes to test for systematic errors. Our final result 
consists of ten different determinations of oip (8.2 GeV), which agree well 
with each other and with our previous results. The most accurate of these, 
when evolved perturbatively to the Z° mass, gives a— (Mz) =.1174 (24). 
We compare our results with those obtained from other recent lattice sim- 
ulations. 

Keywords: Strong Coupling Constant, Lattice QCD, Quarkonium, Per- 
turbation Theory 
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1 Introduction 



Precise measurements of the strong coupling constant a s are important not only 
for strong-interaction phenomenology, but also in the search for new physics. 
Any discrepancy between low-energy and high-energy determinations of this 
coupling could signal the existence of supersymmetry or other phenomena be- 
yond the Standard Model. No significant discrepancies have yet been ob- 
served 0] ; more stringent tests of the Standard Model require further improve- 
ments in precision. In an earlier paper || we showed that lattice simulations 
of quantum chromodynamics (QCD), when combined with the very accurate 
experimental data on the T meson spectrum, provide among the most accu- 
rate and reliable determinations of a s . T's probe the strong interactions at 
the relatively low energies of 500-1000 MeV, where supersymmetry or other 
new physics has little effect. Thus it is important to compare the couplings 
obtained from lattice QCD with those obtained from high-energy accelerator 
experiments, where effects due to a more fundamental underlying theory would 
be much more important. And it is essential that these couplings be measured 
as accurately as possible, with realistic estimates of the uncertainties involved. 
In this paper we review our earlier determination of the coupling, and update it 
to take advantage of new results from third-order perturbation theory, as well as 
new simulations which substantially reduce some of our Monte Carlo errors. We 
also report on several new simulations that further bound our systematic errors, 
particularly with respect to contributions from quark vacuum polarization. 

As discussed in Q , there are two steps in our determination of the coupling 
constant. The first is to create a numerical simulation that accurately mimics 
QCD dynamics. We do this by tuning the bare masses and coupling in a lattice 
QCD simulation until it reproduces experimental results for the orbital and 
radial excitations of T mesons. We use the T family because it is one of the few 
systems for which both accurate simulations and accurate experimental data 
are available. 

Having tuned our simulation, the second step in our determination of the 
coupling is to use the simulation to generate nonperturbative Monte Carlo 
"data" for a variety of short-distance quantities. Comparison with the per- 
turbative expansions for the same quantities then fixes the value of the QCD 
coupling constant. We use the expectation values of small Wilson loops as our 
short-distance quantities. These are very easy to compute in simulations. They 
are also completely euclidean and very ultraviolet, and therefore largely free of 
hadronization or other nonperturbative corrections. Finally, small Wilson loops 
have very convergent perturbative expansions that are known through second 
order for arbitrary rif, the number of light-quark flavors, and through third 
order for rif = 0. 

In this paper we examine each of these steps in detail. We begin in Section 2 
by describing how we tune the simulation parameters. The most important of 
these for our analysis is the bare coupling constant, or equivalently the lattice 
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spacing, used in the lattice QCD Lagrangian. The number and masses of light 
quarks entering through vacuum polarization is also important; we present new 
simulation results that bear on these parameters. In Section 3 we describe 
several different determinations of a-^ using different Wilson loops. Each of 
these sections deals extensively with potential systematic errors. Finally, in 
Section 4 we summarize our results and discuss future directions. 

2 Tuning of the Simulation 

2.1 Procedure 

Given a lattice spacing a, the QCD parameters that determine T properties 
are the bare coupling constant g\ at in the lattice Lagrangian, the bare mass M° 
of the constituent b quarks, and the bare masses wP q of the light quarks that 
enter through quark vacuum polarization. Only the u, d and s quarks are light 
enough to contribute to vacuum polarization appreciably. These parameters 
all vary with the lattice spacing. In a simulation, they must be tuned so that 
physical quantities computed in the simulation agree with the corresponding 
experimental values. The tuning procedure is much simpler, and therefore more 
reliable, if one uses physical quantities that are very sensitive to one of the 
parameters and insensitive to the others. 

Our main interest in this paper is the coupling constant, and so we are 
particularly careful in tuning the bare coupling. We use the mass splittings 
between radial and orbital excitations of the T for this purpose. These splittings 
are ideal since they are almost completely insensitive to the 6-quark mass. The 
spin-averaged mass splittings between IP and IS levels, and 2S and IS levels 
are observed experimentally to vary by only a few percent between the T and 
ip systems, even though b quarks are roughly three times heavier than c quarks. 
This striking insensitivity to the mass of the constituents is an accident, but is 
confirmed by simulations for a range of masses near the b mass. 

These splittings are also quite insensitive to the masses of the light quarks. 
These contribute through vacuum polarization, and affect hadronic masses in 
two ways. First, they allow decays to multi-hadronic final states; mixing with 
these states shifts the masses of the original hadrons. T decay rates are typically 
.1% or less of the mass splittings, and the states we examine are all far below 
the BB threshold. Thus we may ignore such effects in our analysis. The second 
effect of vacuum polarization is to renormalize the gluonic interactions between 
the constituents of the hadron. The typical momentum qr exchanged between 
the b quarks in an T is from .5 to 1 GeV. This is small compared to the c, b and 
t quark masses, and we may ignore their contribution to vacuum polarization. In 
contrast, the u, d and s quarks are effectively almost massless at these energies 
and must be included in a realistic simulation. At the same time, because their 
masses are small relative to qx, our simulation results depend only weakly on 
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their exact values. 

For sufficiently small masses, the dependence of an T mass splitting should 
be linear || : 

AM w AM I 1 + constant X ^ ^ + • • • I , (f ) 

i u,d,s ^ J 

where the renormalized s mass is 50-f00 MeV Q, and the u and c? masses are 
20 or 30 times smaller and therefore negligible. It is very costly to simulate 
lattice QCD with realistic u and d masses. Here that is unnecessary. The 
simple dependence of AM on mP q means that we obtain realistic results if we 
set all three light-quark masses equal to m® s = m°/3, which generates the same 
correction to AM as two massless quarks and a strange quark. Thus m ff = 15- 
30 MeV, and Eq. (|l]) suggests that the dependence on light-quark masses is a 
few percent or less of the total mass splitting, comparable to the Monte Carlo 
statistical errors in our analysis.]] 

There are several other properties of the T system that make it ideal for 
tuning the bare coupling. These mesons are essentially nonrelativistic; the use 
of a nonrelativistic effective action || to exploit this allows a large portion of 
the spectrum to be computed efficiently and precisely Q . They are physically 
small — three or four times smaller than light-quark hadrons — and so do not 
suffer from finite- volume errors even on modestly sized lattices. Finally, we have 
detailed phenomenological quark models that are well-founded theoretically and 
that give us unprecedented control over systematic errors. 

In addition to the bare coupling constant, we must also tune the bare masses 
of the b quark and of the light quarks. We tune the bare 6-quark mass M° by 
requiring that the T mass in the simulation has its correct value of 9.46 GeV. 
Ref. Q presents a detailed discussion. The light-quark masses are tuned until 
the pion and kaon masses are correct. As discussed, we need only the s-quark 
mass, as we set all w/ = 3 light-quark masses to m®/3. 

Finally we note that it is customary in tuning lattice simulations to switch 
the roles of the lattice spacing and the bare coupling constant. Rather than 
choose a lattice spacing and then tune the bare coupling constant to its correct 
value, it is far simpler to choose a value for the bare coupling constant gi a t , and 
then compute the corresponding lattice spacing a using simulation results. All 
explicit dependence on the spacing can be removed from the simulation code by 
expressing dimensionful quantities in units of a or a . The spacing is then not 
needed as an input to the code, but is specified implicitly through the input value 
for gut, or equivalently through (3 = 6/gf at . We determine a from the T mass 

1 It is conceivable that the linear term in Eq. (|l]), which is due to chiral symmetry breaking, 
is strongly suppressed for tiny mesons such as the T, and becomes nonleading. Then the 
dependence on m e fj would be quadratic, with the correct value for m e fl = m,g/\/3. The 
sensitivity to m e fj would then be far smaller and probably negligible for our analysis. 
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splittings AM. The simulation produces these in the dimensionless combination 
a AM; to obtain a, we divide by the experimentally measured values for AM. 

The lattice spacing is a crucial ingredient in our determination of the renor- 
malized coupling a s . As we discuss in Section |3|, the short-distance quantities 
we study specify a s (q*) for a specific value of aq* . The expectation value of the 
axa Wilson loop, for example, gives a s (q*) for q* = 3.4/a. For this to be useful, 
we need to know q*, and therefore a -1 , in physical units such as GeV. Conse- 
quently, the next section focuses on how precisely we are able to determine the 
lattice spacing corresponding to a given value of (3. 



2.2 Results: a 1 Determination 

Our lattice simulations used the standard Wilson action for the gluons, and 
the staggered-quark action and the Hybrid Molecular Dynamics algorithm for 
the light quarks. We employed a nonrelativistic formulation of quark dynamics 
(NRQCD) for the b quarks || 0, M . The n/ = Q gauge-field configurations used in 
our Monte Carlo calculations were provided by G. Kilcup and his collaborators 
(/? = 6,6^1) §, J. Kogut (/3 = 6) [0, and by the UKQCD collaboration (/3 = 
5.7, 6.2) JO]]. The nj — 2 configurations are from the SCRI Lattice Gauge Theory 
Group and their colleagues in the HEMCGC collaboration (/3 = 5.6) jL2|, and 
from the MILC collaboration (f3 = 5.415,5.47)^3). Unfortunately, we were 
unable to obtain configurations with rif = 3 light-quark flavors, which is the 
correct number for T physics. Consequently, we performed complete analyses 
for Uf — and n/ = 2 and extrapolated our results to n/ = 3. The extrapolation 
was the last step of our analysis, and is described in Section |^. 

As discussed above, we use mass splittings in the T system to determine 
the lattice spacing. Specifically, we use two different mass splittings to make 
two independent determinations of the lattice spacing. One is the splitting 
AM(T' — T) between the T' and the T, and the other is the splitting AM(xb — T) 
between the spin average of the Xb mesons and the T. These can be measured 
accurately in a simulation [?j , and are known very accurately from experi- 
ments. Table |l| summarizes the parameters used in our main simulations and 
the results for these two splittings. Our most reliable results are based on the 
/3 = 6 and 5.6 simulations. We use results from the other simulations, including 
those for the splitting between the spin-averaged \c mesons and the spin aver- 
age of the J/t/j and r/ c mesons, to calibrate systematic errors. Our /3 = 6.2 result 



agrees with that in 1 14 



Several factors contribute to the uncertainty in our determination of the 
lattice spacing. We used the lattice NRQCD formalism to simulate heavy- 
quark dynamics ||, and included all relativistic corrections through 0(v 2 ) and 
all finite lattice-spacing corrections through 0(a 2 ). The leading finite-a error 
is due to 0(a 2 ) errors in the gluon dynamics. We estimate this effect using 
perturbation theory M , which indicates that only S states are affected and that 
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am° ff 




splitting 


a AM 


aAM g 


a- 1 (GeV) 


6.0 







1.71 
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-T 


.174(3) 
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.200(12) 


-.005 
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3.15 


Xc - 

Xb 
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.383(10) 
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.0125 


0.80 
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.359(14) 
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-.008 
-.017 


1.30 (5) (20) (5) 
1.44 (6) (4) (6) 
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.05 


0.80 
2.8 


Xc - 

Xb 


-T 


.335(15) 
.307(12) 






6.2 







1.22 


Xb 


- T 
-T 


.124(5) 
.175 (8) 


-.001 
-.0003 


3.58(15)(5)(0) 
3.22(15)(5)(0) 


6.4 







1.00 


Xb 


- T 


.107(16) 


-.002 


4.19 (63) (6) (0) 



Table 1: Simulation results for meson mass splittings a AM and inverse lattice 
spacings a -1 , in GeV, for a range of couplings (i, light-quark masses mj and 
heavy-quark masses M°. The gluonic a 2 corrections aAM g shown are added to 
a AM to obtain the corrected splitting. The error estimates for a -1 are for sta- 
tistical errors, a 2 and i> 4 errors, and errors in the light-quark mass, respectively. 
Experimental values for AM are .563 GeV for T' - T, .440 GeV for X b - T, 
and .458 GeV for \c — ip/Vc- 
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our measured S-state energies should be shifted by 

a AM g = J- (a M q f a AM HFS , (2) 

where AMhfs is the hyperhne spin splitting of the state and M q is the heavy- 
quark mass. We assume 1.5 GeV and 5 GeV for c and b quarks, respectively. The 
corrections we used are listed in Table [l], as are our final values for the inverse 
lattice spacing a -1 . We allow for a systematic error of ±AM g /2 in AM when 
computing the error in a -1 , although our analysis in Q suggests a much smaller 
uncertainty. Note that AAf(T' — T) is almost unaffected by this correction. 
Relativistic corrections of order v 4 are most likely negligible for the T since the 
v 2 corrections, which we include, shift our mass splittings by less than 10%; 
we include a systematic uncertainty of ±1% for this. The v 4 corrections are 
certainly larger for ^'s, where, for example, the J/ip-r] c splitting is a v 2 effect 
and 25% of the Xc—ip/Vc splitting. This suggests v 4 errors could be of order ±6% 
for ip's. Recent simulations |l5| indicate that certain spin-dependent v 4 terms 
can shift levels by as much as 60 MeV, which is 15% of the splitting. We include 
a systematic uncertainty of ±15% for v 4 errors in the tp splitting. 

Our simulations confirm that the 6-quark mass has very little effect on cither 
of the T splittings. The /3 = 6 results show that a 17% change in M° leads to 
changes of only a few percent in the T' — T and Xb — T splittings. (Note 
that the statistical errors in the splittings for different M°'s are correlated. 
Consequently, the statistical errors in the differences between the splittings are 
somewhat smaller than those for any individual splitting.) Since we determine 
M° to within 6% || , the resulting uncertainty in the determination of the lattice 
spacing is probably no more than a percent, which is much smaller than the 
statistical errors. 

Uncertainties in the light-quark mass can also affect our lattice spacing de- 
termination. In our (3 = 5.6 simulations we expect amP s to be somewhere in the 
range 0.01-0.02. This can be inferred from the dependence of the pion mass on 
m q , and allows for uncertainties due to quenching and finite-a errors. Thus we 
want light-quark masses amP q =amP cS in the range .003-. 006. We have simu- 
lation results for ara° fi = .01 and .025. By fitting formula (0) to these results 
we find that am q = .01 should give the correct result to within 4%, which 
equals our statistical error. The correct range of light-quark masses in our 
P — 5.415, 5.47 simulations is roughly amP eS = .005-. 015. We have simulation 
results for a m° ff = .0125 and .050, and again the 6-7% shift caused by changing 
m e ff is roughly the same as our statistical errors for both ip and T splittings.^] 
Note that ip's should be more sensitive to small quark masses than T's since 

2 To compare mass splittings at /3 = 5.415 with those at /3 = 5.47 one needs the expectation 
value of the plaquette at each beta; see the following section. From the plaquette values one 
finds that the lattice spacing at the larger beta is about 12% smaller. Since the a AM's are 
only 5—6% smaller at the larger beta, the splittings AM themselves are actually about 6-7% 
larger for the larger mass. 
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they are roughly twice as large; we saw no evidence for this in our simulations. 
These results all indicate that the m c g dependence is too small compared to 
our statistical errors to allow an accurate measurement .0 This also means that 
the tuning errors associated with a m e ff are no larger than our statistical errors, 
and so we take our statistical errors as a measure of the uncertainty due to this 
parameter. 

We checked for finite- volume errors by computing the charmonium splittings 
using lattices that are 1.5 fm and 3.0 fm per side. We observed no difference, 
indicating that these errors are smaller than the 2% statistical errors in these 
tests. The lattices we used at p = 6 and 5.6 are both 16et« 1.35 fm per side; the 
T's are half the size of the ip's, with a radius of about .2 fm. We therefore expect 
finite volume errors in our mass splittings that are substantially less than 2%. 

We estimated the electromagnetic shifts of the T masses using a potential 
model. For individual mesons, we found mass shifts of approximately 1 MeV, 
with smaller shifts for the splittings between them. These are too small to affect 
our result. 

Our final values for a _1 's are listed in Table [j], obtained by dividing the ex- 
perimental values for the splittings AM by the corrected Monte Carlo simulation 
results a AM + aAM g . The error estimates for the a _1, s include statistical er- 
rors in a AM, as well as systematic errors associated with the finite-a correction 
aAMg, v A corrections, and the light-quark mass mP q . Other systematic errors 
are negligible. 

Perhaps the most striking feature of these simulation results is the significant 
disagreement at /3 = 6 between a -1 computed using the X' — T splitting and that 
computed using the \b — T splitting. Taking proper account of correlations, this 
disagreement is five standard deviations: our simulation gives 1.43(3) for the 
ratio of these splittings, rather than the experimental value of 1.28. Thus the P = 
6 simulation is inconsistent with experiment. This is because in this simulation, 
in contrast to nature, n/ = 0; there is no light-quark vacuum polarization. The 
disagreement is much smaller when n/ =2, as is apparent in the y3 = 5.6 data. 
And, as we will demonstrate, it disappears completely when we extrapolate rif 
to three. 

As expected, using an incorrect value for rif leads to inconsistencies such as 
the one found in our /3 = 6 simulation. Perturbation theory, though not justified 
at the momenta relevant for these systems, provides a qualitative explanation for 
this discrepancy. The centrifugal barrier makes the average separation between 
the quarks in the P state Xb larger than for the S state T or T', as is familiar 
from hydrogen or positronium. As a result, the typical exchanged momentum 
for Xb quarks, q Xb , is smaller than gx'- The perturbative binding energy is 
given by a 2 s (q)Cp Mb/16, with q = qr> for T' and q Xb for Xb- Since q Xb < qr>, 
the Xb is more tightly bound. However, for rif = 0, this effect is exaggerated, 

3 This insensitivity to m e g is because m c g is so small in our simulations. Our rif = 
simulations are equivalent to m ff = oo and give results that are quite different from rif = 2. 
So shifts would become apparent, even with our statistics, for sufficiently large m e g . 
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Figure 1 : NRQCD simulation results for the spectrum of the T system, including 
radial excitations. Dashed lines indicate experimental values for the triplet S 
states and for the spin average of the triplet P states. The energy zero from 
simulation results is adjusted to give the correct mass to the T(l 3 Si). Results 
are from a simulation with rif — (filled circles) and from one with rif = 2 
(open circles), using a -1 = 2.4 GeV for both. The errors shown are statistical; 
systematic errors are of order 20 MeV or less. 



as a ( ° } (q) increases more quickly than af 1 (q) with decreasing q. Thus, for 
rif < 3, AM(xb — T) should be underestimated relative to AM(T' — T), as is 
observed. Fitting to data would then require a larger a -1 for AM(\b — T) than 
for AM(T' - T). 

We end this section by displaying in Figures |] and [2] results from the /3 = 6 
and 5.6 simulations for several of the low- lying excitations and spin splittings, 
compared with experimental values. The agreement is excellent and supports 
the reliability of our simulations. We emphasize that these are calculations 
from first principles; our approximations can be systematically improved. The 
only inputs are the Lagrangians describing gluons and quarks, and the only 
parameters are the bare coupling constant and quark masses. In particular, 
these simulations are not based on a phenomenological quark potential model. 
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Figure 2: NRQCD simulation results for the spin structure of the lowest-lying 
P states. Dashed lines indicate experimental values for the triplet P states. 
Masses are relative to the spin-averaged state. Results are from a simulation 
with rif = (filled circles) and from one with nj = 2 (open circles), using 
a~ 1 = 2.4 GeV for both. The errors shown are statistical; systematic errors are 
within about 5 MeV. 

3 Determination of the Renormalized Coupling 

3.1 The Coupling Constant from Wilson Loops 

Having tuned the simulation, we performed Monte Carlo simulations to generate 
"data" for a variety of short-distance quantities. We determined the coupling by 
matching the perturbative expansions for these quantities to the nonperturba- 
tive Monte Carlo results. For short- distance quantities we chose the expectation 
values W mj „ of Wilson loop operators. In the continuum, 



where P denotes path ordering, is the QCD vector potential, and the inte- 
gral is over a closed ma x na rectangular path. Loop operators for small paths 
are among the most ultraviolet, and therefore most perturbative, objects that 
can be studied in lattice QCD simulations. Unlike most other quantities used 
to determine the QCD coupling, the loop operators are truly short-distance 
quantities in euclidean space. There are no corrections for hadronization, and 
nonperturbative effects are expected to be very small. For example, the leading 
nonperturbative contribution to W m ,n due to condensates is probably from the 




(3) 
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gluon condensate, with 



6W m , n = -*-^pl { a s F>). (4) 

Most studies find that (a s F 2 ) is of order .042 GeV 4 Q. Since a -1 ranges 
from 1.2 to 4.2 GeV in our simulations, we expect condensate contributions to 
— lnWi,i, for example, to be in the range of .2-. 01%, much too small to be 
important here. When rif ^0 there are also contributions from quark conden- 
sates, but these are suppressed by a 2 s and so are probably even smaller. The 
tiny size of such effects make the W m , n for small m and n ideal quantities for 
determining the coupling in lattice QCD, particularly given the ease with which 
they can be computed in simulations. 

To obtain four independent determinations of the coupling, we used expec- 
tation values for the four smallest loops on the lattice: the plaquette Wi } 2, 
W%,3, and W2,2- Each of these loop operators is very different from the others; 
as different, for example, as various moments of a structure function. Each 
is affected differently by nonperturbative effects and higher-order uncalculated 
perturbative corrections. The contribution of the gluon condensate, for exam- 
ple, is 16 times larger for than for Wi t %. By comparing results obtained 
from different loop operators we can bound such systematic errors. 

Each of our expectation values has a perturbative expansion of the form 

-laW&tf =^2^ lf \m,n)(a^\q m , n )y, (5) 

i=l 

where a P is a new nonperturbative definition for the coupling constant intro- 
duced in our earlier paper Q to facilitate lattice calculations. The scale q m ,n is 
the average gluon momentum in the first-order contribution to W m>n , computed 



directly from the Feynman diagrams as described in 1 17, |18| . 

In Table |^ we list the perturbative coefficients through third order for rif = 0, 
and through second order for rif = 2 |l9| . Unfortunately, the n / dependence of 
the third-order coefficients has not yet been computed. Given that the second- 
order coefficients depend only weakly on rif by design jD], it is likely that 
the rif = third-order coefficients are also good approximations when rif = 2. 
We assume this in our analysis, but when estimating errors at rif — 2 we take 
the size of the entire rif = third-order contribution as an estimate of the 
uncertainty due to rif dependence. When rif = 0, we estimate the truncation 
error in perturbation theory to be of order a P (q m ,n) times the leading order 
contribution. 

Note that the plaquette W\ i has no third-order contribution. This is because 
the coupling a P is defined in terms of the plaquette ; the absence of third and 
higher-order corrections is merely a consequence of our conventions. Truncation 
errors in the plaquette's expansion reappear when our coupling is converted to 
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loop 


Cl 


C2 




C3 


t* y m , n 






n f = 


fly = 2 


ny=0 




-lnWi,i 


4.19 


-4.98 


-5.57 





3.40 


-lnWi^2 


7.22 


-7.57 


-8.51 


2.6 


3.07 


-lnWi^ 


10.07 


-9.60 


-10.89 


5.3 


3.01 


- lnW 2 .2 


11.47 


-10.58 


-11.84 


11.1 


2.65 



Table 2: Coefficients for the perturbative expansions, in powers of a P {q m ^ n ), of 
small Wilson loops. Scale q m ,n is the average momentum carried by the gluon 
in the first-order correction. 

more standard couplings, such as Om§: 

agv)(Q) = a£'\ ( &*Q){l + 2c£' ) /* (6) 

Here the third-order coefficient X-^ w 0.95 for ny = 0j|^]. The third-order 
coefficient is new since our first paper. Unfortunately, the ny dependence of this 
coefficient is not known. However, the variation of this coefficient as ny goes 
to two or three is unlikely to be large. The factor e 5 / 6 in the scale is chosen 
to eliminate n / dependence in the second-order coefficient of the expansion Jl8[ | , 
and therefore also removes much of the ny dependence in third order. As above, 
we use the ny = value for Xjjg throughout our analysis, but when ny = 2 we 
take the size of the entire third-order term as our estimate of the uncertainty 
due to ny dependence. 

The coupling a P was defined to coincide through second order with the 
continuum coupling ay defined in Jl8| , [HJ from the static-quark potential. 
The third-order correction to the static-quark potential has recently been com- 
puted |2lJ , leading to 

ap\Q)=a^\Q) {l + X v ( ) 2 + • • •} , (7) 

where Xy — 1-86 — .14 ny + X^, which is 2.81 for ny = 0. Note that this 
expansion has infrared divergences in fourth-order and beyond, due to residual 
retardation effects in the static quark potential j2^] . 

3.2 Results: a P Determinations 

Monte Carlo simulation results for the expectation values of the Wilson loop op- 
erators are summarized in Table || ^3|. We also tabulate the values of a P {q mjl ) 
obtained by matching perturbation theory to Monte Carlo simulation results. 
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.1584(6) 
1 6^7 r«S 


5.6 


2 


.025 
.010 
.010 
.010 
.010 


-lnWi.i 
-InWi.i 

-lnPt^ 3 
- In T¥ 2) 2 


.5719(0) 
.5709 (0) 
1.0522(1) 
1.5123(2) 
1.8337(3) 


.1792(0) 
.1788 (0) 
.1828 (30) 
.1832 (40) 
.1907(80) 


5.7 







-In 


.5995(0) 


.1829(0) 


5.415 


2 


.0125 


-lnWi,i 


.6294(0) 


.2075 (0) 


5.47 


2 


.050 


-In Wi,i 


.6134(0) 


.1993 (0) 


6.2 







-InWi.i 


.4884(0) 


.1398(0) 


6.4 







-lnWi,i 


.4610(0) 


.1302 (0) 



Table 3: Expectation values of Wilson loop operators for small loops, and the 
corresponding a P 's for a variety of lattice QCD parameters. The uncertainties 
listed for the expectation values are Monte Carlo statistical errors. Those listed 
for the ap's are estimates of the truncation errors in perturbation theory. 



The uncertainties quoted are our estimates of the potential truncation errors in 



perturbation theory; see Section 3.1. The only other potential sources of error 
are nonperturbative effects, and as discussed, these are almost certainly negligi- 
ble compared to truncation errors. Finite-volume errors are much less than 1% 
for such small loops. 

The values for the various coupling constants in this table are all different. 
This is because the coupling-constant scales q m , n are different for each operator 
and for each parameter set. To compare these results we must first evolve 
the running coupling constants to a common scale. In Table we present the 
couplings evolved to 8.2 GeV, which is the scale we chose in g]. To generate 
these values, we converted the corresponding q m ,nS from units of a -1 to GeV 
using the lattice spacings inferred from each of the T or ip mass splittings for 
which we have simulation results. We then evolved the couplings to 8.2 GeV by 
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loop 




ai n/) (8.2GeV) 












Xb - T 




Xc ~ ip/Vc 


6.0 







-InWi.i 

1 TT7 

In M/1,2 
-lnWi, 3 
- In Vp2,2 


.1552(10)(0) 
.1556 (10)(6) 
.1560(H)(6) 
.1565(H)(8) 


.1505(H)(0) 
.1509 (H)(6) 
.1512(H)(6) 
.1517(H)(8) 




5.7 







-InWi.i 


.1528(18)(0) 




.1465 (61) (0) 


6.2 







-InWi.i 


.1569 (23) (0) 


.1519(23)(0) 




6.4 







-InWi.i 


.1515(67)(0) 






5.6 


2 


.010 
.010 
.010 
.010 


-lnW M 
-1nWi i3 
- In ^2,2 


.1794 (24) (0) 
.1777(24)(30) 
.1770(24)(40) 
.1767(23)(71) 


.1781(33)(0) 
.1764 (32) (30) 
.1757 (32) (40) 
.1754(32)(71) 




5.415 


2 


.0125 


-lnWi,i 


.1748 (34) (0) 




.1696 (78)(0) 



Tabic 4: Values of a P (8.2GeV) from several operators W m ,n and a variety of 
tunings for QCD simulations, with different (3's, n/'s, and meson mass splittings 
used to fix a -1 . The two uncertainties listed are due to uncertainties in the 
inverse lattice spacing, and to truncation errors in the extraction of a P using 
perturbation theory. 

numerically integrating the evolution equation for ct P . We used the universal 
second-order beta function together with the n/ = third-order term for a P . 
The n / dependence of the third-order beta function is unknown, but the entire 
third-order term generally has negligible effect. This is especially true for our 
most important results at (3 = 6 and 5.6, since 8.2 GeV was chosen to be very 
close to the q m ^ n 's and very little evolution is required. 

If one groups the various couplings in this table according to the splitting 
used to tune the simulation and the number of light-quark flavors n / , one finds 
that the values within a single group arc completely consistent. In particular, 
results obtained using different loops are in excellent agreement, which shows 
that our estimates of the errors caused by truncating perturbation theory are 
reasonable. Also, the coupling constants obtained from the plaquette using /3's 
ranging from 5.7 to 6.4, corresponding to scales q\ t \ ranging from 4.8 GeV to 
14.2 GeV, agree well. This demonstrates that the evolution of our coupling con- 
stant dp is well described by the perturbative beta function; no lattice artifacts 



14 
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(GeV) 



Figure 3: Values of the QCD coupling constant ap determined from the pla- 
quettc in simulations with differing lattice spacings corresponding to (3 = 5.7, 6, 

6.2 and 6.4, all with n/=0. The coupling constant is plotted versus the average 
momentum q\ t \ carried by gluons in the plaquette at the various lattice spac- 
ings, with qi } i = 3.4/a. The line shows the coupling constant evolution predicted 
by third-order perturbation theory. 

are apparent. This is also illustrated by Figure |^, where we plot the coupling 
constant otp{qi,i), obtained from the plaquette, versus the effective momentum 
scale (7^1 = 3.4/0 at which the coupling is measured on each lattice. The simu- 
lation results for the running of a P agree well with the prediction of third-order 
perturbation theory . 

3.3 Extrapolation to rif = 3 

The coupling constants in Table [| from simulations with different n^'s are sig- 
nificantly different, as are the couplings from simulations tuned using different 
meson mass splittings. Our final step is to extrapolate to rif =3, which is the 
correct number of light-quark flavors for T and tp physics. The extrapolated 
results, which are shown in Table should all agree, and do. To make the 
extrapolation, we paired rif = and nf = 2 simulations as indicated in the table. 
For each separate combination of Wilson loop and meson mass splitting, we 
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(3 loop ag' (8.2 GeV) 

Xb ~ T T' - T Xc - VV»?c 



6.0,5.6 



2.2 



1,3 



1,1 



1.2 



1946 (41) (0) 
1913(42)(52) 
1897(42)(69) 
1889(40)(120) 



1960 (61)(0) 
1927(58)(54) 
1912(58)(71) 
1903(57)(125) 



5.7,5.415 



InWi 



1884 (57) (0) 



■ 1841(146)(0) 



Table 5: Values of ap 3 ' (8.2 GeV) from different operators and different tunings 
of the QCD simulation. The two uncertainties listed are due to uncertainties 
in the inverse lattice spacing, and to truncation errors in the extraction of a P 
using perturbation theory. 

extrapolated l/a P using the corresponding q p 's from the two simulations. 

We chose to extrapolate l/a P rather than a P because numerical experiments 
using third-order perturbation theory suggest that l/a P is significantly more 
linear in n/. To see how the couplings from our simulations might depend on 
rif, note that the T splittings that we use to determine the lattice spacing probe 
QCD at momentum scales qr of the order 0.5-1 GeV. Thus when we choose a 
lattice spacing that gives these splittings their correct physical values, we are in 
effect tuning the QCD coupling constant in our simulation to have its correct 
value at the scale qr- (If rt/^3, the simulation's coupling will have the correct 
value only at q~r.) This means that the couplings in our n/ =0 and 2 simulations 
agree with the correct n/ = 3 coupling at qy. 



This equation specifies the dependence of the couplings obtained in our sim- 
ulations onnj, but we are unable to use it directly since perturbation theory 
is not particularly reliable at q-f- Nevertheless, we can use this relation to test 
different schemes to extrapolate rif as follows. Taking qy = 1 GeV, we set all 
the couplings at that scale equal to some large value, say .65. We then evolve 
all three to 8.2 GeV using the three-loop beta function. Finally, we compare the 
8.2 GeV coupling extrapolated from nj — and 2 with the nj — 3 coupling ob- 
tained by evolving from qr . Extrapolating a P gives results that are "correct" to 
within 1.4%, while extrapolating l/a P is correct to within 0.3%. This exercise 
indicates that we should extrapolate the inverse coupling and that the extrap- 
olation errors are probably less than 1%. Such errors are negligible relative to 
the other systematic and statistical errors. Nevertheless, it would be desirable 
to repeat our analysis using simulations with n/ = 3 or even n/ = 4. 



OL 



x ?(qy) = a?(q r ) = a?(q r ). 



(8) 
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Eq. (||) played a key role in the earliest determinations of the running cou- 
pling constant using lattice QCD p!| . These studies used only nj — simula- 
tions. As can be seen from our results, the coupling at n/ =0 is 25% smaller 
than the correct n/ = 3 coupling. This correction was estimated in these earlier 
papers by perturbatively evolving the n/ = coupling down to qr , changing n f 
to three, and then evolving back up to the original large scale, which is 8.2 GeV 
in the present analysis. This procedure suggests a correction of 15-20%, which 
our simulations show to be an underestimate but within the error range quoted 
in the earlier papers. We emphasize that there is no inconsistency between 
these earlier analyses and ours. Our simulations with rif = give results that 
are identical with the earlier work. What is different here is that we have actual 
simulation results at n/ ^ and so get to nf = 3 using extrapolation, rather 
than a perturbative analysis that is well-motivated but only partly justified. 
That the sizable correction due to light-quark vacuum polarization was so ac- 
curately predicted using perturbation theory strengthens our confidence that 
our nonperturbative treatment of vacuum polarization is correct. Note that if 
we use the perturbative analysis to correct just our n/ = 2 couplings, ignoring 
our n/ = couplings, we obtain results that are in excellent agreement with the 
extrapolated coupling |^6| . 

Our final results for a P in Table | agree well with each other and with 
our earlier results Q. In particular, the 5<r discrepancy between results using 
different T splittings at nf — Q disappears completely at n/ =3. This is highly 
nontrivial; we are in effect counting the number of light-quark flavors that affect 
real upsilons. It provides confirmation that the quark vacuum polarization is 
correctly included in our simulations and extrapolation. 

3.4 Conversion to a^s 

To compare with nonlattice determinations of the coupling constant, we have 
converted our results to the MS definition of the coupling, using Eq. (||) with 
Ams = .95 ± -95. Our results are listed in Tabic |g| an d together with our a P 's in 
Table |, are the main result of this paper. The MS results are somewhat larger 
than in our earlier paper because we now use the n/ = value for X-^ , rather 
than setting it to zero as before. Our estimate in the earlier paper for the size 
of this term was correct and was included as an error. Consequently, our old 
results are consistent with our new results within errors. 

To further facilitate comparisons with other analyses, we have numerically 
integrated the third-order perturbative evolution equation for a-^ and applied 
appropriate matching conditions at quark thresholds [^7| to evolve it to the mass 
of the Z°. The results for our ten determinations are shown in Table 0. For 
matching we assumed MS masses of 1.3(3) GeV and 4.1(1) GeV for the c and b 
quarks respectively ||, |27| . The uncertainties in these masses can shift the final 
coupling constant by less than half a percent; we ignore them. 
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loop g<|(3.56 GeV) 

Xb ~ T T' — T Xc - V"/»?c 



6.0,5.6 


-InWi.i 


.2258(56)(74) 


.2277(83)(75) 




-In Wi, 2 


.2213(56)(99) 


.2232 (79) (102) 




-InWi's 


.2192(57)(116) 


.2212(79)(119) 




- In W-2,2 


.2181 (54)(176) 


.2200 (77) (183) 


5.7,5.415 


-lnW hl 


.2174(76)(67) 


.2117(197)(62) 



Table 6: Values of o£^(3.56GeV) from different operators and different tunings 
of the QCD simulation. The two uncertainties listed are due to uncertainties in 
the inverse lattice spacing, and to truncation errors in the extraction of a P and 
conversion to a^s using perturbation theory. 






loop 












Xb - T 




Xc - ip/Vc 


6.0,5.6 


-lnWi,i 
-lnWi' 2 
-lnWi^ 3 
- In W 2 ,2 


.1174(15)(19) 
.1163(15)(26) 
.1157(15)(31) 
.1154(14)(46) 


.1180(22)(20) 
.1168(21)(27) 
.1162(21)(31) 
.1159(20)(48) 




5.7,5.415 


-lnWi.i 


.1152(20)(18) 




.1136(52)(16) 



Table 7: Values of a^(Mz) from several operators and various tunings of the 
QCD simulation. The two uncertainties listed are due to uncertainties in the 
inverse lattice spacing, and to truncation errors in perturbative expansions. 
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4 Discussion and Conclusions 



In this paper we have demonstrated that lattice simulations provide among 
the simplest, most accurate, and most reliable determinations of the strong 
coupling constant. Our ten different results, tabulated in Tables ^-0, are in 
excellent agreement with each other. Indeed, all but one of them agree with 
our best determination to within its uncertainty; that is, to within the smallest 
error bars. Our best result implies 

0.3706 (288) for Q = 1.3 GeV w M c and n f = 3 
0.3701 (288) for Q = 1.3 GeV « M c and n f =4 



0.2234 (93) for Q = 4.1 GeV « M b and n f = 4 (9) 
0.2233 (93) for Q = 4.1 GeV w M b and n/ = 5 



I 0.1174 (24) for Q = 91.2 GeV = M z and n/ = 5 , 

with errors due to lattice-spacing and perturbation-theory uncertainties com- 
bined in quadrature. These results are about 1 er higher than our previous 
results [|J . The shift is entirely due to the new third-order term in the perturba- 
tive formula, Eq. (^|), relating the lattice coupling a P to a-^. Our Monte Carlo 
simulation results are essentially identical to those in our earlier paper. The 
shift relative to our earlier result is only 1 a because we previously estimated 
the size of this third-order term accurately. 

The bulk of our effort in this analysis was devoted to understanding and 
estimating the systematic errors. We varied every parameter in the simulation. 
We used four different short-distance quantities to extract the coupling, and 
three different (infrared) meson splittings, in two different meson families, to 
tune the bare coupling or lattice spacing. We demonstrated that the gross 
features of T and ip physics are accurately described by our simulations. We 
explored the role of light-quark vacuum polarization for a range of light-quark 
masses. Our simulations were sufficiently accurate to show that n/ = is the 
wrong number of light-quark flavors for T's. Only when we extrapolated to 
nf = 3, the correct value, did our simulation results agree with experiment. To 
see how robust our results are, we redid the analysis but with various ingredients 
missing. The corresponding shifts in a^(Mz) are listed in Table |§ omitting 
the n/ extrapolation led to the only appreciable difference. 

The various parts of our analysis agree well with the results of other groups. 
The a P 's that we extract from Wilson loop operators agree to within statistical 
and truncation errors with those obtained by very different techniques [29[| . This 
is the easy part of the analysis. The remainder, involving the determination of 
lattice spacings, has now also been duplicated. A recent analysis of simulation 
results from the Fermilab and SCRI groups, both of which employ a totally 
different formalism for 6-quark dynamics, gives a^(Mz) = -116 (3), in complete 
agreement with our results ]3C| . 
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omit 0{a 2 ) gluonic corrections 
omit tadpole improvement of NRQCD 
omit 0(v 2 ,a,a 2 ) corrections in NRQCD 
omit extrapolation (use n/ = 2) 



-0.6% 
-0.5% 
+0.9% 
-4.7% 



Table 8: Changes in the coupling constant at Mz when different parts of our 
simulation or analysis are omitted. 

Lattice coupling constant determinations such as ours enjoy a fundamen- 
tal advantage over traditional methods based on perturbative high-energy pro- 
cesses, allowing significantly greater accuracy. The systematic uncertainties in 
the perturbative parts of the analyses are similar in both approaches, but the 
nonperturbative elements differ substantially. When we tune our simulation to 
reproduce the T spectrum, we are in effect directly tuning the QCD scale pa- 
rameter A-j^r. Consequently, a 5% simulation error in a mass splitting results in 
a 5% error in A^g, which implies only a 1% error in a^(Mz)- In high-energy 
determinations, however, one measures the coupling constant rather than the 
scale parameter, and usually only through small radiative corrections to an 
clcctroweak process. Measuring A^ is intrinsically much more accurate than 
measuring Ojjg. 

There are prospects for substantially improving the accuracy of our result 
fairly soon. We list sources of error in our value for a^L{Mz) in Table ||. 
The dominant error is due to truncation in perturbative expansions, specifically 
because the n / dependent parts of our third-order coefficients have not yet been 
calculated. The agreement we observe between couplings from different loop 
operators, each with its own perturbative series, suggests that our estimates of 
this systematic error are realistic or even pessimistic. Nevertheless, our total 
error could be cut in half by computing this nt dependence, particularly for 
Eq. (||). This is a straightforward perturbative calculation. For this paper, 
we halved our statistical errors for our n/ = simulations; the same should 
be done for rif ^ 0. Use of an improved gluon action would remove the need 
for the a 2 correction in the Xb ~ Y analysis, while it already has negligible 
effect on T' — T. The additional cost would be small Using a relativistic 
formulation of c-quark dynamics, rather than NRQCD, might allow accurate 
results from the charmonium spectrum. A simulation with either rif = 3 or 4 
light quarks would eliminate the extrapolation error and would require perhaps 
only twice the computational effort needed for rif — 2. Finally, simulations with 
larger light-quark masses m e ff would allow us to pin down more accurately the 
dependence on this parameter. 

Our lattice determinations of the strong coupling constant agree well with 
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Uncertainty 



Unknown rif dependence in third-order perturbation theory 1.9% 

Statistical error in determination of a -1 .9% 

Light-quark masses .9% 

Extrapolation in n/ .3% 

Finite a and 0{v A ) errors .2% 

Fourth-order evolution of a-^ .01% 



Table 9: Sources of error in our best determination of oi^JMz) 



most determinations based on perturbative high-energy processes. This fact 
provides striking evidence that the nonperturbative QCD of hadronic confine- 
ment and the perturbative QCD of high-energy jets are the same theory. 
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